Cargar las librerías

library(tidyverse)
library(sf)
library(mapview)
library(GADMTools)
library(ggspatial)
library(leaflet)
library(leaflet.extras2)
library(spdep)
library(spatstat)  
library(raster)
library(smacpod)
library(ggspatial)
library(DT)

Caso de estudio

Cargar los datos

covid <- read_csv(url("https://zenodo.org/record/4915889/files/covid19data.csv?download=1"))
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   FECHA_CORTE = col_date(format = ""),
##   UUID = col_character(),
##   DEPARTAMENTO = col_character(),
##   PROVINCIA = col_character(),
##   DISTRITO = col_character(),
##   METODODX = col_character(),
##   EDAD = col_double(),
##   SEXO = col_character(),
##   FECHA_RESULTADO = col_date(format = ""),
##   rango_edad = col_character(),
##   lon = col_double(),
##   lat = col_double()
## )
covid_p <- covid %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)
covid_p %>%
  filter(FECHA_RESULTADO=="2020-12-11") %>%
  ggplot() +
  geom_sf()

m_p<-covid_p %>% filter(FECHA_RESULTADO == "2020-12-10") %>% 
  mapview(layer.name="puntos")

m_p
peru <- gadm_sf_loadCountries("PER", level=3)
lima_sf <- peru %>%
  pluck("sf") %>%
  # Filtramos los datos espaciales solo de Lima metropolitana
  filter(NAME_2 == "Lima") %>%
  # Editamos algunos errores en nuestros datos espaciales
  mutate(NAME_3 = ifelse(NAME_3 == "Magdalena Vieja",
                         "Pueblo Libre", NAME_3))
mexico<-gadm_sf_loadCountries("MEX",level=2)
cdmx_sf<-mexico %>%
  pluck("sf")
mexico
## $basename
## [1] "./"
## 
## $sf
## Simple feature collection with 1854 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -118.3689 ymin: 14.53292 xmax: -86.71014 ymax: 32.71804
## Geodetic CRS:  WGS 84
## First 10 features:
##      ISO NAME_0          NAME_1              NAME_2    TYPE_2    ENGTYPE_2
## 3683 MEX Mexico  Aguascalientes      Aguascalientes Municipio Municipality
## 3665 MEX Mexico  Aguascalientes            Asientos Municipio Municipality
## 3687 MEX Mexico  Aguascalientes            Calvillo Municipio Municipality
## 3652 MEX Mexico  Aguascalientes               Cosío Municipio Municipality
## 3690 MEX Mexico  Aguascalientes         Jesús María Municipio Municipality
## 3675 MEX Mexico  Aguascalientes Pabellón de Arteaga Municipio Municipality
## 3657 MEX Mexico  Aguascalientes     Rincón de Romos Municipio Municipality
## 3661 MEX Mexico  Aguascalientes  San José de Gracia Municipio Municipality
## 3663 MEX Mexico  Aguascalientes            Tepezalá Municipio Municipality
## 3336 MEX Mexico Baja California            Ensenada Municipio Municipality
##                            geometry
## 3683 MULTIPOLYGON (((-102.448 21...
## 3665 MULTIPOLYGON (((-102.229 22...
## 3687 MULTIPOLYGON (((-102.613 21...
## 3652 MULTIPOLYGON (((-102.2537 2...
## 3690 MULTIPOLYGON (((-102.5781 2...
## 3675 MULTIPOLYGON (((-102.4122 2...
## 3657 MULTIPOLYGON (((-102.19 22....
## 3661 MULTIPOLYGON (((-102.4122 2...
## 3663 MULTIPOLYGON (((-102.2376 2...
## 3336 MULTIPOLYGON (((-114.7158 2...
## 
## $level
## [1] 2
## 
## $hasBGND
## [1] FALSE
## 
## attr(,"class")
## [1] "gadm_sf"
covid_count_mx <- covid %>%
  group_by(DISTRITO, FECHA_RESULTADO) %>%
  summarise(casos = n()) %>%
  ungroup() %>%
  complete(FECHA_RESULTADO = seq.Date(min(FECHA_RESULTADO, na.rm =T),            max(FECHA_RESULTADO, na.rm = T),                                      by="day"),
           nesting(DISTRITO), fill = list(n = 0))
## `summarise()` has grouped output by 'DISTRITO'. You can override using the `.groups` argument.
covid_count <- covid %>%
  group_by(DISTRITO, FECHA_RESULTADO) %>%
  summarise(casos = n()) %>%
  ungroup() %>%
  complete(FECHA_RESULTADO = seq.Date(min(FECHA_RESULTADO, na.rm =T),
                                      max(FECHA_RESULTADO, na.rm = T),
                                      by="day"),
           nesting(DISTRITO), fill = list(n = 0))
## `summarise()` has grouped output by 'DISTRITO'. You can override using the `.groups` argument.
covid_sf <- lima_sf %>%
  mutate(DISTRITO = toupper(NAME_3)) %>%
  full_join(covid_count, by = "DISTRITO", "FECHA_RESULTADO")

class(covid_sf)
## [1] "sf"         "data.frame"
covid_sf %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf()

m_sf <- covid_sf %>% 
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  mapview(layer.name = "distritos")

m_sf

Múltiples capas

ggplot() +
  geom_sf(data = covid_sf %>% 
            filter(FECHA_RESULTADO == "2020-12-11")) + 
  geom_sf(data = covid_p %>% 
            filter(FECHA_RESULTADO == "2020-12-11"))

m_p +m_sf

Visualización de datos espaciales

Patrones de puntos

covid_p %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf(aes(col = SEXO), alpha = .2) +
  facet_wrap(.~SEXO)

p1<-covid_p %>% 
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  mapview(layer.name = "points", zcol = "SEXO", burst = T)
p1
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:raster':
## 
##     select
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Dos o más variables

covid_p %>%
  filter(FECHA_RESULTADO == "2020-04-11" |
         FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf(aes(col = SEXO), alpha = .2) +
  facet_grid(SEXO~FECHA_RESULTADO) +
  guides(col = F)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

m1 <- covid_p %>%
  filter(FECHA_RESULTADO == "2020-04-11") %>%
  mapview(zcol = "SEXO", layer.name = "2020-04-11 - SEXO")

m2 <- covid_p %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  mapview(zcol = "SEXO", layer.name = "2020-12-11 - SEXO")
m1 | m2

Composición

covid_p %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf(data = covid_sf) +
  geom_sf(aes(col = EDAD), alpha = .2) +
  scale_color_viridis_c(option = "B") +
  annotation_scale() +
  annotation_north_arrow(location = "tr",
                         style = north_arrow_nautical)+
  theme_bw()

Datos en polígonos

Una variable

covid_sf %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf(aes(fill = casos))

covid_sf %>% 
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  mapview(layer.name = "casos", zcol = "casos")

Dos o más variable

covid_sf %>%
  filter(FECHA_RESULTADO == "2020-04-11" |
         FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf(aes(fill = casos)) +
  facet_grid(.~FECHA_RESULTADO)

d1 <- covid_sf %>%
  filter(FECHA_RESULTADO == "2020-04-11") %>%
  mapview(zcol = "casos", layer.name = "2020-04-11 - casos")

d2 <- covid_sf %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  mapview(zcol = "casos", layer.name = "2020-12-11 - casos")
d1 | d2

Composición

covid_sf %>%
  filter(FECHA_RESULTADO == "2020-12-11") %>%
  ggplot() +
  geom_sf(aes(fill = casos)) +
  scale_fill_viridis_c(option = "F", direction = -1) +
  annotation_scale() +
  annotation_north_arrow(location = "tr",
                         style = north_arrow_nautical)+
  theme_void()

Variación espacial del riesgo

covid_subset <- covid %>%
  filter(FECHA_RESULTADO == "2020-05-05")

covid_win <- owin(xrange = range(covid_subset$lon),
                  yrange = range(covid_subset$lat))
covid_ppp  <-  ppp(covid_subset$lon, 
                   covid_subset$lat, 
                   window = covid_win)
densidad_raster_cov <- raster(density(covid_ppp, bw.ppl), 
                              crs = 4326) %>%
  mask(lima_sf)
densidad_raster_cov %>% 
  mapview()
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

Detección de clústers

Datos puntuales

Estadísticas de escaneo espacial (Spatial Scan Statistics-SSS)

covid_subset_posi <- covid %>%
  filter(FECHA_RESULTADO == "2020-05-05") %>%
  mutate(positividad = ifelse(METODODX == "PCR", 1, 0))

covid_scan_ppp <- ppp(covid_subset_posi$lon, 
                      covid_subset_posi$lat,
                      range(covid_subset_posi$lon),
                      range(covid_subset_posi$lat),
                      marks = as.factor(covid_subset_posi$positividad))
covid_scan_test <- spscan.test(covid_scan_ppp,
                               nsim = 49, case = 2, 
                               maxd=.15, alpha = 0.05)
covid_scan_test
## $clusters
## $clusters[[1]]
## $clusters[[1]]$locids
##  [1] 1103 1371 1092 1375 1381 1104 1380 1378 1091 1074 1097 1376 1089 1076 1373
## [16] 1379 1096 1087 1088 1095 2093 2122 2117 2102 2106 1105 2121 1107  605 2086
## [31] 1098  592 1101 2111  596 2096 1070  599 2120  597  971 2094 2119  607
## 
## $clusters[[1]]$coords
##          [,1]      [,2]
## [1,] -77.0472 -12.11374
## 
## $clusters[[1]]$r
## [1] 0.02709301
## 
## $clusters[[1]]$pop
## [1] 44
## 
## $clusters[[1]]$cases
## 110344 
##     31 
## 
## $clusters[[1]]$expected
## [1] 11.17098
## 
## $clusters[[1]]$smr
##   110344 
## 2.775046 
## 
## $clusters[[1]]$rr
##   110344 
## 2.873837 
## 
## $clusters[[1]]$propcases
##    110344 
## 0.7045455 
## 
## $clusters[[1]]$loglikrat
## [1] 20.05829
## 
## $clusters[[1]]$pvalue
## [1] 0.02
## 
## 
## 
## $ppp
## Marked planar point pattern: 2316 points
## Multitype, with levels = 0, 1 
## window: rectangle = [-77.18195, -76.64829] x [-12.468713, -11.634232] units
## 
## attr(,"class")
## [1] "spscan"
# Construimos el centroide del clúster
cent <- tibble(lon = covid_scan_test$clusters[[1]]$coords[,1],
               lat = covid_scan_test$clusters[[1]]$coords[,2]) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = F)  

# Construimos el área del clúster en base al radio
clust <- cent %>%
  st_buffer(dist = covid_scan_test$clusters[[1]]$r)
cluster <- mapview(clust, alpha.regions = 0, color = "red") 

points <- covid_subset_posi %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  mapview(zcol = "positividad", alpha.regions = .4, alpha = 0) 

cluster + points 

Datos agregados (en polígonos)

Autocprrelación espacial (global):Moran

covid_sf_subset <- covid_sf %>%
  filter(FECHA_RESULTADO == "2020-05-05") %>%
  mutate(casos = replace_na(casos, 0))
covid.nb <- poly2nb(covid_sf_subset, queen=TRUE,snap = 0.13)
covid.lw <- nb2listw(covid.nb, style="W", zero.policy=TRUE)
moran.test(covid_sf_subset$casos, covid.lw)
## 
##  Moran I test under randomisation
## 
## data:  covid_sf_subset$casos  
## weights: covid.lw    
## 
## Moran I statistic standard deviate = 3.0548, p-value = 0.001126
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.090491242      -0.023809524       0.001399995

Autocorrelación espacial (local):Getis Ord

breaks <- c(-Inf, -1.96, 1.96, Inf)
labels <- c("Cold spot",
            "Not significant",
            "Hot spot")
covid_lg <- localG(covid_sf_subset$casos, covid.lw)

covid_sf_lisa<-covid_sf_subset %>% 
  mutate(cluster_lg=cut(covid_lg, include.lowest = TRUE,
                        breaks = breaks, 
                        labels = labels))
covid_sf_lisa %>%
          ggplot() + 
          geom_sf(aes(fill=cluster_lg)) +
          scale_fill_brewer(name="Clúster", 
                            palette = "RdBu", direction=-1) +
  theme_bw()